Skip to main content
Internet Archive's 25th Anniversary Logo

Full text of "DIRECT CALCULATION OF SUBLIMATION ABLATION"

See other formats


> 



NASA TECHNICAL TRANSLATION 



NASA TT F-l4,529 



DIRECT CALCULATION OF SUBLIMATION ABLATION 



H. Kubota 



Translation of: "Shoka abureshon no 
chokusetsu kalho", Koku Uchu Gijutsu 
Kenkyusho Hokoku, [Technical Report 
of National Aerospace Laboratory] , To K I O^ 

No. 239, June, 1971, PP- 1-25 



SlSftfon ""fee? «u.. ^"^ ^° " "fo. 



G3/33 



N72-32936 



Unclas 
U1014 



;♦' 



fjiSSzi?^ 



kO> 



>S J^nZ'S?: 




NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
WASHINGTON, D.C. 205^6 AUGUST 1972 



Reproduced by' 



SPnngf,e|JvA%i?,'«» 



DIRECT CALCULATION OF SUBLIMATING ABLATION 
Hlrotoshi Kubota 

Abstract. A direct and rapid calculatlonal . 
method of sublimating ablation for a hypersonic '■ 
heat-shielding device is presented. This is 
based on a numerical analysis of the Hartree- 
Womersley method which reduces parabolic partial 
differential equation to a corresponding 
ordinary one by means of a kind of implicit. 

1. Introduction 

When a spacecraft reenters the atmosphere, powerful shock waves /L*** 
are produced around the craft, and the air currents behind the 
shock|waves assume an extremely high temperature. It has been made 
clear in recent research that an effective way to protect the craft 
from this sort of aerodynamic heating is to adopt mass transfer 
cooling, in which a relatively small quantity of cooling gas or 
cooling liquid is blown off into the external fluids. 

This utilizes chiefly the following three effects: (1) By mix- 
ing the low-temperature cooling gas with the main currents of high- 
temperature gas, the air-current temperature near the surface of the 



^Received March 23, 1971- 

an 

First Aerodynamics Department, 

***Numbers in the margin ^indicate pagination in the original foreign 
textf. 



object Is lowered. (2) The temperature gradient Is reduced on ac- 
count of the Increased thickness of the boundary layer due to the 
Increased mass Inside the boundary layer, and also due to variations 
of the velocity distribution and temperature distribution Inside the 
boundary layer. (3) Heat absorption occurs because of chemical re- 
actions or phase changes Inside the boundary layer or on the surface 
of the substance. 

Possible heat-protection methods are transpiration cooling, film 
cooling, local mass Injection cooling, or ablation cooling. They 
each have their advantages and disadvantages, and It Is difficult to 
determine which of these methods would be most effective under all 
possible conditions. Nevertheless, the ablation cooling method con- 
tains the above-mentioned effects of (2) and (3), and there Is no 
doubt that It still remains one of the most promising heat-protection 
modes as long as one Ignores the actual, economical problems such 
as repairs of the eroded surface of the object. It should be noted 
that this method has a rather high cooling effect. 

This method was based on the idea of protecting at all costs the 
most Important parts of the craft, while allowing, or rather posi- 
tively encouraging, the melting of the surface substances of the 
craft by means of external heat. Research concerning this was per- 
formed at the early stages by Sutton [1] and Roberts [2]. During 
the process of ablation, generally speaking, a liquid phase boundary 
layer and a gaseous phase boundary layer are formed on the surface 
of the object. That is, the solid phase surface of the object melts /2_ 
producing a liquid-phase boundary layer, over which a gaseous bound- 
ary layer is formed consisting of a mixture of the evaporated gas 
and the external gas. The thermal energy in the high-temperature 
air currents is transmitted to the surface of the object through 
these boundary layers. Part of it is absorbed into the latent heat 
of evaporation and the latent heat of melting, while the remainder^ — 
is absorbed in the form of the thermal capacity of the structural 
materials. = 



The substances chiefly used as the surface materials are high 
polymer plastics, graphite, etc. Other materials are also known 
In which protection from heat Is attained by bringing about a chemi- 
cal reaction on the surface or by creating a char layer Inside the 
object [3, 4, 5, 6, ?]• 

If the object Itself Is subjected to sublimation, the surface 
substance will change directly from a solid Into a gas . Therefore, 
there will be only a gaseous boundary layer on the surface of the 
object. Nevertheless, this does not Interfere essentially with the 
validity of the analytical method, which Is based on a two-phase 
boundary layer Including also a liquid boundary layer. Thus, this 
paper will deal with sublimating ablation and will clarify Its 
mechanism. 

The writer and his colleagues already on a previous occasion 
used the micro-disturbance theory to analyze sublimating ablation of 
hemispherical objects, using Teflon materials [5]. The results 
showed that the calculated values differed considerably from the 
test values when there was a relatively low stagnation point tem- 
perature. It would be significant to resort to an electronic com- 
puter and use a unified method of obtaining a direct, numerical 
solution, not only in order to check out the results obtained pre- 
viously, but also to give the theory a more universal validity. From 
this standpoint, the questions dealt with in this report may be sum- 
marized in terms of the following three points: 

(1) Essentially, the ablation phenomenon ought to be determined 
by the process of chemical reactions of the surface substances in 
combination with the aerodynamic heat transfer process in the circum- 
ference of the object. Thus, if the conditions inside the substance 
and the conditions of the main current are given, the physical vari- 
ables on the surface of the substance, particularly the ablation 

speed, the surface temperature, and the heat transfer coefficient, 
will be determined categorically. 



(2) The flow becomes dissimilar downstream from the stagnation 
point of axially symmetrical objects. The; physical variables in 
the flowwise direction are sought by applying the Hartree-Womersley 
method, proposed by Hartree and Womersley [8] and further developed 
by Smith et al. [9] for the dissimilar terms in the boundary layer 
equation describing this area. The purpose of this method (here- 
after abbreviated .as| the H-W method) is to reduce a parabolic partial 
differential equation to a corresponding ordinary one by means of a 
kind of Implicit scheme .[ 

The least square error method is applied to the two-point bound- 
ary value problem of the ordinary differential equation obtained by 
the H-W method as described above, and a general solution is proposed. 
A simplified version of this method has been used the Nachthelm and 
Swigert [10] in cases when there are as many as two unknown boundary 
conditions. However, in the present case, seven departure values 
are necessary to carry out numerical integration of the ordinary 
differential equation. Of these, only one Is given explicitly, and 
the other six are unknown. These six departure values (unknown 
boundary conditions) are related by three bo ns train t| conditions. 
The remaining necessary number of boundary conditions are given at 
the Itnternal edge|of the boundary layer. Two methods are used in 
problems of this type: the initial value method, and the quasi- 
linearization method. In the initial value method [11], the departure 
values (unknown boundary conditions) are suitably assumed, the 
numerical values are integrated, and iterative calculations are con- 
tinued until the conditions of the external edgel of the boundary 
layer are satisfied. In the quasi-linearization method [12], the 
original equation is linearized, and a non-homogeneous solution — 
consisting of a lower-order approximate solution — is sought. Both 
methods give good results when the number of departure values (unknown 
boundary conditions) to be assumed is small, but there is a tremen- 
dous increase in the difficulty when the number becomes larger Rur- 

thermore the presence of iconstraintf conditions on the wall also 
complicates the problem. 



The method used here is a further development of the initial 
value method. No matter how the number of unknowns or the number 
of restraining conditions ;are| increased, the method can be applied 
in principle as long as there is agreement with the fundamental 
equation. It is also believed that the method can be applied to 
other problems as well, as long as they are boundary layer type 
problems (with asymptotic behavior at the upper limits of the domain) , 

Let us assume a semi-stationary ablation, in which the surface 
sweepback velocity, and thence the ablation velocity, is almost un- 
affected by the time. (For a detailed discussion, refer to [5]-) 
The thermal properties and physical quantities of most ablation sub- 
stances are not clear [7]. Teflon (polytetraf luoroethylene) is one 
of those substances which are subject to sublimating ablation (it 
has been observed experimentally that there is probably a melting 
process, but since the temperature range of this process is extremely 
small], it is assumed that there is only a sublimation process) and 
for which the physical quantities have been clarified. Thus, Teflon 
is used as the model in these calculations. 

In the method described here, as long as the state of the main 
current and the state inside the object are given, it is possible 
to carry out automatically the above-mentioned processes (1), (2), 
and (3), and the ablation physical quantities can be calculated in a 
unified manner as far as any point downstream from the stagnation 
point of a blunt-headed object. Such methods of direct solution 
usually rely on the memory capacity of a computer, but in the H-W 
method the upstream solutions are discarded as soon as they are used. 
This makes operations less difficult, and this is one of the advan- 
tages of the method. 



2. Notation 

a: velocity of sound 

C: Chapman- iRunf sln| number 13- 

Cp: specific heat at constant pressure 

Cv: specific heat at constant volume 

D: two-component diffusion coefficient 

d: boundary layer thickness 

E: square error In two-point boundary value problems 

H: total enthalpy function 

H „„: effective ablation heat 
ef f 

K: concentration function 

L: latent heat of sublimation 

M: molecular weight 

M : Mach number of main stream ! 

cx> ' 

m: ablation velocity 

Pr : Prandtl number 

p: pressure . - - - 

q: heat transfer rate 

R: gas constant 

Rb : radius of curvature of object 

Re: Reynolds number 

r^ : radius of curvature of cylinder 

Sc : Schmidt number 

(s, n ) : transformed coordinate system 

T: temperature 

t: time 

(Uj V): velocity components In spatially fixed system of coordinates 

(u, v) : velocity components In system of coordinates fixed to object 

V, : surface sweepback velocity due to ablation 

(X.j.y): spatially fixed system of coordinates ^ ; 

(xj y): system of coordinates fixed to object 'i 

X.: boundary conditions 

z: transformed concentration function 



specific heat ratio of air 

error in two-point boundary value problem 

criterion of square error 

heat transfer coefficient 

cooling effect due to ablation 

coefficient of viscosity 

density 

transformed i g^tream functionl 

stream function] 

transformed enthalpy function 



Suffixes: 



a 
b 
e 

, g 
n 

no 
s 

sh 
t 
w 

1 



( )' 
( ~) 

_ K 

( ) 



state with ablation 

state inside the object 

state at external edgejof boundary layer 

state of the gas phase 

n^^ state 

state' with no ablation 

state of solid phase 

state immediately after shock waves 

state at stagnation point of main strea: 

state on wall of object 

state at stagnation point of object 

state of ablation gas 

state of external air 



state of main stream] 

differential 

dimensional quantity 

quantity for non-dimensionalization 



Basic Equations 



Let us consider semi-stationary ablation, in which the surface 
"sweepbackl velocity is more or less uniform and does not depend on 
the time. It is assumed that all variations caused by ablation are 



J7 



confined within the boundary layer and have)no effects on the exter- 
nal current. There are no chemical reactions Inside the boundary 
layer, and the boundary layer consists of a two-component mixed gas 
(the gas caused by ablation and the external air). The thickness of 
the boundary layer is assumed to be sufflclentr y small in comparison! 
with the radius of curvature of the object. 



In a flow field such as that shown in Figure 1, the compressible 
laminar boundary layer equation for a hemispherical object (Appendix 
A) can be written as below after 
transformation of the independent 
variables by means of the Less- 
Dorodnitsyn transformation: ;; 



Shock waves 



Shock wave, layer] 



5= I pelitu.rMx, Ji=p,u,r,{2s) * I £-dy 
Jo J op, 

'(C5S")'+/3(-e.'-(5>'+l)2}+Vi"(5> + ,)=i?i, 
+ Wu,' |C(1-A)5>"(?>' + 1)}' =i?2, 



^&*'^'+*'('*''"'')=^'- 



(3 


1) 


(3 


2) 


(3 


.3) 


(3 


4) 



Te{X),Ue(X) 

'Boundary 
layer 



. M., n^ 




Figure 1. System of coor- 
dinates and flow pattern 



/4 



Here, ( ) ' is the differential re- 
lated to n, and 



a=.l-CNCM^^^ 



(3.5) 



C, Pr, and Sc refer to the Chapman-Rube sin number, the Prandtl num- 
ber, and the Schmidt number (Appendix B) ; they are each functions 
of (|) , ^3 and z. Here, the values of cj)', ^, and z are the following 



f=u/u.-l, V^H/H.-l z=K 



(3.6) 



The stream function fi is defined as follows:] 



3» _„„, 9a 



fl=V 2s {^is,,,)+,,^ 



(3.7) 



Therefore, cj) is expressed as follows: 

=i=r' ^ / (3.8) 



If we use an equation of state for a perfect gas, the following 
will apply because of the condition that there are no pressure vari- 
ations in the y direction in the boundary layer: 

f-={i- ^^-^>^}4r.| (3.9) 

- 

Due to the equation for the total enthalpy (A.l6), the temperature 
will be 

T _ !r+l-|«.W(^'+l)' l , 



In the basic equations (3-2) - (3.^), the dissimilar terms 
R-, , Rp, and R^ on the right side will each be: 






(3.11) 



j?,=2.{(^'+ i)^i^-r'ii|±£)} J ■( 3 : 1 2 ) 



(3.13) 



The dissimilar terms will be expressed by these items and by the 
differential term of s. The treatment of this will present one of 
the problems. That is, because of the presence of dissimilar terms 
R-, , Rp, and R^, it will be impossible to express the flow field by 



means ;of o ne s imilarity! parameter n , and a separate solution must 
be sought for each station (that is, for the various values of s). 

4. Boundary Conditions 

Since the basic equations (3.2) - (3.4) are, respectively, 
three-step, two-step, and two-step differential 'equations with^Te'- 
spect to ()) , ^i*, and z, altogether seven boundary conditions are neces- 
sary. However, only four conditions are explicitly given: the con- 
dition that there is no slipping on the wall (u, = 0), and the 

W 



)' 



condition that there is connection with the external flow at the 
end of the boundary layer (u ^ u^, H •> H^, z -> ; y ^ y^)- When 

these are written in transformed form, they are as follows: 

,= ; 5i'=-i. I (4.1) 

^ = oo ; 5>'=0, 1^ = 0, 2=0 I (4.2) 

However, the following thre]e Icons traint i conditions occur on 
account of the physical conditions on the surface of the object where 
ablation takes place. 

(i) Relationship between, wall temperature T^ 
and mass loss rate m 

When an ablation substance^ is_ considered, __lts_mass-_loss__. rate 
\m is a characteristic function determined by the wall tempera- 



ture T and the pressure p. It can be written as follows: 
w 

m = FCTu,.J)\ (4.3) 

However, various test results and data [13] have shown that in the 
low-temperature region the effects of pressure are so small as to 
be almost negligible. Rashis and Hopko [l4] have obtained this 
relationship for Teflon as shown in Figure 2. When this is approxi- 75. 
mated in a quadratic equation, it becomes: 



m = — ptVi= (pv)w 



(4.4) 



One can thus obtain an approximate equation to take the place of 
equation (4.3) • 

On the other hand, from the aerodynamic viewpoint, if the state 
of the main ■streamj is given, the relat ionship between the wall tem- 
perature T and the mass loss m can be described by the]- .. 



following equation: 

^=G(T,o, 0)1 (4.5) 



10 



7.0 



S56.0 



E 
V5.0 



^4.0 



3.0 



2.0 



1.0 



0.0 



« * 



5.0 



Si 4.0 



E 
(J 



,3.0 



X 

IS 



500 



550 600 



650 



700 7501 

ru<'c) 



2.0 



1.0 



/?6= 0.75cm 




i? t= 2.00cm 
3.00 



Figure 2.. Chemical rela- 
tionship between the wall 
t emp e ra ture an d the ma s s 
loss rate for Teflonl 
[TiT] . 



500 



550 600 



650 700 
Tw{'C) 



750 



Figure 3. Aerodynamic rela- 
tionship between the wall 
temperature and the ablat-lon 
rate for Teflqnj 



Here, Q Is a quantity which Is 

determined according to the Mach 

number of the main ,s treaml , the state of the. stagnatJ:On point of the 



T andm 
w 



are ' unknowns 



main stream), the state of the object, etc 

In Equation (4.4) and Equation (4.5), and It Is anticipated that the 
two unknown numbers will be determined uniquely with respect to 
various ambient conditions by converting these two equations Into 
a simultaneous equation. In actual fact, the relationship in Equa- 
tion (4.5) was obtained as in Figure 3 according to [5]- By over- 
lapping Figure 2 and Figure 3 diagrammatlcally , it is possible to 

obtain the solution with T and m at the respective intersection. 

w 

In this report, we proceed one step further and attempt to—at- — 
tain a comprehensive solution, in which Equation (4.4) is also in- 
cluded as part of the boundary conditions. When Equation (4.4) 



C 



11 



is written In transformed fashion. It assumes the following 
appearance: 



X {Ai 4- AzT.T ,0 + Ja (J I r,o)2J 



(4.6) 



(11) Mass equilibrium on the surface of 
the object 



K. , the concentration of the various component gases on the 



iw 



surface of the object, cannot be designated without relation to the 
-Wall velocity v . When v is determined, k. will also be deter- 

f| WW IW 



mined. Generally speaking, v^ , the velocity of the 1 gas in the y 

D dK,l 

direction, is the sum of the diffusion velocity "a", ITjand the con- 



vection velocity v. Therefore, 



pK,v, = -pD^J^+pKiV 
dy 



(4.7) 



Only the gas components produced by sublimation will be blown awayj 
from the surface' of the object. Therefore, one can posit the physi- 
cal condition that there is no movement of the external air Inwards 
or outwards along this surface. Under this condition, the v^^ with 

respect to the external air (suffix 2) will be 0. If we use the 



fact' that Kp = 1 - K, we obtain the following: 



Pvo 



H-f)- 



H-p. V. d-K", )=0 



(4.8) 



Equation (4.8) can be rewritten as follows in transformed 
fashion: 



m. r^=^-2.(^), 



(4.9) 



(1). 



Ata constant stagnation point temperature, if there is an increase 



of T there will be a decline of the temperature gradient inside the 
w 

boundary layer. Consequently, the thermal input towards the surface 

of the object will decline, and there will be a smaller ablation 

velocity. 



12 



(ill) Energy equilibrium on the surface 
of the object 

On the surface of the object where ablation Is taklijg place. 



part of the heat quantity V dyJz 



which Is about to be transferred 



through the boundary layer to the object will be absorbed by the 
latent heat of ablation, and the remainder of It will be transferred 
Inside the object. Thus, If we take Into account the conduction of 
heat In the x direction Inside the object, the energy equilibrium on 
the surface of the object can be written as follows: 



('41. =^^>~>-+-{(^);.. 



dy /». 






(4.10) 



The second Item on the right side Is obtained as the solution of the 
heat conduction equation for the Interior of the object: 



_,3'T fT\ dT_ _:._ n^ 



(4.11) 



Nevertheless, since Equation (4.11) is an elliptic partial differen- 
tial equation. It would be Impossible to combine It In a simultaneous 
equation with a parabolic boundary • layer equation and apply the H-W 
method. Furthermore, with the range of stagnation point temperatures 
being handled here, heat diffusion In the x direction Is negligible 
In comparison with that In the y direction [5]- Therefore, Equation 
(4.11) can be written as follows: 

dT 






dy 



(4.12) 



Here, m = - p^ v^ = ( p v)^ 



/6 



Let us solve Equation (4.12) under the following boundary 
conditions : 

7=0 ; f = T„ -1 / 

y=-yr.T^Tjj (4.13) 

In dimensionless form, the following solution will be obtained: 

r=c, exp(^^)+J (4.14) 



13 



Here , 



c,= (r.-r.)/{exp(— -|L^.)-i|^ 

Cz=T„—ci, 



R 



VS. 



(4.15) 



Therefore, the following will apply 

dT\ _ P 



(3T\ P l\-l. 

(4f)., .=P' /""•''•' (4f)., . 



(4.16) 



By substituting Equation (4.16) in Equation (4.11), we will obtain 
the following equation for energy equilibrium in transformed fashion: 



1 — aZm ' PalCtoUtTeTt 



Here 



{-^ gg-/+W^-+M 



(4.17) 



E\w 



/=Z/C,2'5",. 



(4.18) 



The wall temperatures T appearing in the various equations can 



be expressed as follows, in accordance with Equation (3.10): 



r. = (!P-. + l)/(l-a2.) 



(4.19) 



They are related to ^ and z . 

^ WW 



5. Calculating Methods 



5.1. Hartree-Womersley Method 



The H-W method, which was described in the Introduction, trans- 
forms into the difference system the dissimilar items R, , R„ and R 

on the right side of the basic equation, which is a parabolic partial 



differential equation, as well as the differential terms] having the 



14 



form of dF/ds which are contained In the boundary conditions. The 
equation Is thus reduced to an Implicit ordinary differential equa- 
tion. The two-point differences of equal Intervals used here were: 

OF _ F-F.-i l (5-1) 

ds As / 

However, the precision will Improve further If three-point differences 
with variable Intervals are used. Nevertheless, the conve^lon velo- 
city distribution on the surface Is smooth In ablation problems. 
Therefore, there Is no necessity for discussing the difference sys- 
tem with reference to the discontinuous convection velocity which 
was dealt,, with In [15], and It will be sufficient to use two-point 
differences . 

According to the calculating diagram In Figure 4, the solution 

for the n point (the suffix n Is omitted In order to make the equa- 
tions easier to read) will Include 

the solutions for points upstream u n-2 n—1 

from It. That Is, at the nth 



point, one will use F„_-| i the known 

solution for the points upstream, 
and the ordinary differential equa- 
tion for n will be solved in the n 
direction. The dissimilar terms 
on the right side of the basic 
equation will be written as follows : 



Solu-I 
tlon 



Fn-ii. 1)) 



Solu- 
tion 

Fn-iiv) 



Solu-] 
tlon 



I 



■^S 



Figure 4. Calculation dia- 
gram for the H-W method 



i?i=-|^{(i' + l)(*'-0'.-l)-«>"(!*-?>.-l)} 

As mentioned above, 4 -,,(!)' -,5*1' -, , z -, are known solutions 

' ^n-1' ^ n-1' n-1' n-1 



(5.2) 



'The stagnation point of a blunt object becomes the point 



of departure for calculations, and there are no solutions for points 
upstream from it. However, with respect to the finite difference 
amplitude at this point. As, it can be stated that 2s/As = 0,. and 
therefore one can obtain a solution corresponding to cases where 
Rj = R2 = R3 = (similar solution) . 



15 



5.2. Two-Point Boundary Value Problem 

The ordinary differential equation obtained above does not have 
the necessary number of boundary conditions on the wall, and It 
therefore becomes ja two-point boundary value problem of the type 
described in the Introduction. 

Let us posit the following in order to give the problem a uni- 
versal validity. 

"In a boundary layer type ordinary differential equation with /7 
n steps (asymptotic behavior), m (m < n) boundary conditions are 
given at the .lower limit of the domain, and (n - m) of the conditions 
are unknown. These (n - m) unknown boundary conditions are linked 
by means of I {I < n - m) constraint] conditions. At the upper limit 
of the domain, (n - m - I) boundary conditions are given. We seek 
a solution satisfying these conditions." 

Let us express the boundary conditions In terms of Xj (j = 1, 
2, ..., n) . First of all, the conditions which are explicitly given 
are X-,, Xp, .... X (their number is m) . In the initial value method, 

(n - m - I) boundary conditions, excluding this and the I constralntj ' 
conditions, must be supposed as the departure values for numerical 
integration. These boundary conditions which are to be supposed as 
expressed as X^^^-j^, \+2' "' "^n-l ^^^^^^ number is] (n - m - Z)]. If 
these (n - m - I) conditions are suppos ed, the remaining I conditions 
will be determined from the ^c ons traint conditions, and therefore the 



solution insidje the domain will be governed by the Independently 

assumed! X , -, , X ^„, ...., X ,. If the solution is expressed as S., 
- - — I m+1 m+il n-i J 

it can be written as follows : 

Sj =Si (Xm-ny x„+2, , x.-d \ (-5.3) 

(y=i. 2, , «)1 



16 



In the present problem, n = 7j m = 1, and I 
which are given explicitly are: 

and the conditions which are to be assumed] are 

^2=5>"^, Xi=W,, Xi=z,\ 



3; the conditions 
(5.5) 



The conditions which are determined via the ;co-nstralnt| condi- 
tions are: 



Jf6=/. =X6(Jt2, ^3, X4) = ( 4 • 9 ) ^ 



(5.6)/ 



The solutions inside the boundary layer are: 

Si=55=Si (X2, X3, X^) 
S2=55' = S2 (X2, Xz, Xt) 

S3^p" = S3(X2, X3, X'4) 

St^'F = St (X2, X3. Xt) 

S6=r=S5(X2, X3, Xt) 

S6=2=S6 (ATa, A'3, Xt) 
S^^^' = St (_X2, X3, Xi) 



(5.7) 



As was'' mentioned before, the number of conditions necessary at the 
end of the boundary layer is (n - m - Z). In the present case, this 
number would be 3. They are given by equation (4.2). However, as 
was mentioned in the Introduction, in these boundary layer type 
equations, all the physical quantities have an asymptotic behavior 
at n ■> «>. Therefore, the following also ought to apply supple- 
mentarlly: 



l? = oo ; ?>" = 0, if' = 0, i'=0 



(5.8) 



If the conditions in Equation (5.8) are not present, there may be 
solutions in which ^a certain quantity will overshoot along the way 
(Figure 5), so that the solution is no longer an asymptotic on e. 



In order to avoid this phenomenon, the supplementary conditions 
in Equation (5.8) become necessary as full conditions. 



17 



Let us use suffixes I and II 
to indicate the first approxima- 
tion and second approximation of 
solutions at the external edge of| 
the boundary layer (the upper 
limit of the domain) by Iterative 

calc ulations . ^jln this case, the 

modified quantities of the de- 
parture values (the supposed un- 
known boundary conditions and the 
boundary conditions determined by 
the ' constraint! conditions) and 
the modified quantities of the 
solutions will be: 




Figure 5' Explanation of 
solutions, with overshoot 

in the center 



^X,= (,XOii-(iXi)i, (.-=2, 3, 4) 
^S/=(Sy )ii -(Sy ) I , (; =2,3, .7) 



Furthermore, according to the yairlatlonal principle 



ASy=2]^AXf (;=2, 3, , 7) 



(5.9) 
(5.10) 



(5.11) 



Here, the operator 



D 



DX. 

1 



is defined as follows : 



^ ydXk d 
DXr'dXi'^LdXi dXk 

According to Equation (5.10) and Equation (5.11) » 

4 

'(Sy)lI=(Sy+2]^^X)l (/=2. 3. 7) 



(5.12) 



(5.13) 



Since Equation (^.2) and Equation (5.8) are the desirable conditions 



at the external edge of the boundary layer of (S.)jj, If <S . is the 
error between (S.)jj and them, it will be: 



'-Sj =Sy +J]^AXi U = 2. 3, , 7) 



(5.1^) 



The suffix II was omit'ted here in order to give universality to the 
approximations. It is clear that the necessary and sufficient condi- 
tion for convergence of the solution is to give 6j a value of 0. 
For this purpose, let us take a square error of 



/8 



18 



L^' 



(5.15) 



and assume that the convergent solution will be when this assumes a 
value of 0. In actual practice, let us take a small positive number 
and adopt the following as the criterion for convergence: 



£ < £ 



(5.16) 



The conditions for the least square error are obtained from 
the following: 



dE 



d (AXi) 



= 0=2, 3, 4) 



In the present problem. 



dE 



1 = 2 



1=2 
4 



g-^=0J:Oe„ + ^e.<AX<=0. 



(5.17) 



(5.18) 



Here 






y=2 

7 



»-i:(S>' . 



(1 = 2, 3, 4, *=2, 3, 4) 



(5.19) 



The following is obtained from the above simultaneous 
equation: 



19 



AA'2= - 



AX3= - 



AA'4= - 



DEM = 



©21 


023 


024 


031 


033 


034 


041 


043 


044 


©22 


021 


024 


032 


031 


034 


042 


041 


044 


022 


023 


021 


032 


033 


031 


042 


043 


041 


011 


012 


013 


021. 


022 


023 


031 


032 


033 



/DEM , 



/DEM 



/DEM , 



(5.20) 



In one trial, a modifier term of AX^ is added to X^, giving X^ + AX^ , 
and we move on to the next trial. 



However, in order to obtain AX. in the manner of Equation (5.20) 

it Is necessary to find, by some method or other, the value in Equa- 
tion (5.19) taking the form of DS./DX. . This is necessary at all the 

stations in the boundary layer. Differentiating the basic equation 
by X. (1 =2, 3, ..., 7), we obtain the disturbance differential 



equation : 



/ s 



^+^A{-^-('''-^»^&^'> 



+ !»' 



„S^_dR\ 






3 f C ,\.3z',., s ,,, # _dK3\ 



(«=2, 3, 



•. 7) 



=0 (i^j) 
(«, y=2, 3, 



•. 7) 



(5.21) 



(5.22) 
(5.23) 



Solving this equation under the following boundary conditions : 



(5.24) 



let us substitute this solution, together with the values of 8X^/8X^ 



(1=2, 3, 4, k=5, 6,7) (known on account of the constraint 



20 



conditions between the boundary conditions). In Equation (5-12) , the 
defining equation of p-^:— . The S. In Equation (5.2^1) Is a function 

determined as In Equation (5-7). All of the boundary conditions 
(departure values )' for the disturbance differential equation obtained 
here, are given as In Equation (5.24). The solution Is obtained by 
making a simultaneous equation of the total of 21 differential equa- 
tions. Including the 3(n - m) disturbance differential equations 
(there are l8 of them In the present problem) and the three original 
equations. The square error Is calculated at each of the calculating 
points (at each of the scale Intervals of the numerical Integration 
In the n direction). Nachthelm and Swlgert [10], established In ad- 
vance a suitable n a. (point for stopping calculations) and Investl- 

StOp -^ Ir-Jr- o 

gated the effects of Its value on the solution. In this report, we 
posit the following criterion: 

|Sy(,)|^a,y| ( 0) . = constant) (5.25) 

J 

Solutions departing from these conditions are regarded as divergent 

solutions, and the calculations are stopped at that point. This 

point is called n = n ^ • Then we proceed to the next trial. The 
^ stop 

end point of the calculating region given in advance is also called 

n = n . . In this way, the convergent solution is one which matches 
' 'stop '^ ' ° 

the conditions of Equation (5.I6) while also fulfilling the condi- 
tions of Equation (5.25). Thus, the n = n satisfying the conditions 

of Equation (5. 16) will correspond to the external edge| of the bound- /g 
ary layer, and it is possible to find also the thickness of the 
boundary layer, which was one of the unknown numbers. 

The flow chart for the calculations is as shown in Figure 6. 



21 




x = x+/ix C 111 ,'") 
s = s +Js 



:« 



;! 



( »• = 2 , 3 , 4 ) 



>7=0 



3?= J7+ JJ? 



(Runge-Kutta-Gill) 



El^ 



_('x. 



\3^ 



J-X-f ( 1=2,3,4) 



V=ve 



^:3Z7 




Figure 6 



1 — beginning; 2 



Flow chart for the H-W method and the method of solution 
for two point boundary value problems 

- read In conditions of main .■jstreamr and Inside the 



assume Initial values of X 



3 — stagnation point; 
- Xs, Xs, and X? are determined from lconstralnt|, 



2.1 J 



X 



3 5 



object ; 

X:^; 5 — Xs, Xs, and X? are determined from lconstralnt|, conditions 
6 —original differential equation and disturbance different lal_ equa- 
tion are solved as a simultaneous equation; 7 — - E Is calculated; 8 — 
E Is written out; 9 — end; 10 — does the solution converge?; 11 — 



(downstream region); 12 — obtained solution is included in Ri , R2 , R 
13 — solution is written out; l4 — AXj_ (1 = 2, 3, 4) is calculated." 



3 5' 



22 



Numerical Calculation Results and 
Considerations of Them 



The Teflon (polytetraf luoroethylene) used in the numerical cal- 
culations is a fluoride polymer, but it decomposes into a monomer 
upon pyrolysis. According to Madorsky [l6], more than 90^ of it is 
CpF. . Therefore, let us here regard C2Fij as the ablation gas and 

assume that the boundary layer consists of a two-component mixed gas 
(air and CpFh gas). The physical values of these substances are as 
follows : 

C2F4 gas: ^^^^00 j 

Cpi=0.32 cal. gr-i. ."K-i/ 

ffi=5. OOA r s 

j2,(2.2)=o.90 / 

(See Appendix B) 



/lO 



Air: 



A/2=29 

Cm=0. 28 cal . gr-i . °K-i 

72=3. 62 A 

.O2«-2)=0.76 



(See Appendix B) 



Solids Teflon : c>.=o.22 cai . gr-i. "k-* 

pi =2. 19 gr .cm-a 

Kt = 6. OOX IO-*''cal.cin->.sec-». °K-> 

L=35kcal • mol-i 



The following are assumed as the conditions of the main stream : 



A/oo=5.74, r,=800°C~1500°C, r.=20''C, 
^( = 0. 5atm~2atm, ^»=0. 75cm~3. QOcni 



In Figures 7(a) and 7(ti) is shown the state of convergence of 
the departure values of a two-point boundary value problem when— the— 
calculations were performed according to the method described in the 
preceding chapter. The calculations were carried out on the HITAC- 
5020 computer at the National Aerospace Laboratory. Figure 7(a) 



shows the state of convergence between (j)'' , 

W 



m 



w 



^w' 



the values 



/ll 



^ 



23 



-1.0 



¥m 



-0.8 



-0.5 

z'w 



--0.4 



|l-0.6 



ii 1 



^0.4 



1 — r 



1 — I — I — r 



, -0.2 



-0.1 



1.0 



'w, 



Zw 



: Departure value s j 
^^'=0.0, 7'»'=0.90,arH-=0.'01 



:rlals 


E 





5.41191 


1 


1.55245 


2 


2.16962 


3 


0.11206 


4 


0.00029 


5 


0.00008 




2.0 



p; 



1.5 



0.5 







1.0 



0.5 



^y 3 4 5 6^ 7 

'No. of trial si 



Figure ?• State of convergence of 
solution : 

(a) convergence of initial value: 

a; = 0, T, = 1000°C, ^. = 1.00cm, ^, = latm. 



tO.5 



-1.0 




Figure ?• State of converg- 
ence of solution: 

(b) velocity profile, 

i*=0, T, = 1060°C, ^ = 1.00cm, i>;=latm.| 



assigned at random at the beginning (number of trials = 0), and 
d) , 1" , z' , determined according to them. The square error E in 

each trial is also shown here. The square error was calculated at 
n = n • As for the divergent solutions departing from the criteria 

of Equation (5-25), they are those at the maximum n position (ffT ) 

S G Op 

as long as the conditions of Equation (5.25) are satisfied. The 
following were adopted as the conditions of Equation (5.25): 



24 



032 


= 1.5 


(for 


*') 


0)3 


= 1.0 


(for 


*") 


0)i| 


= 0.5 


(for 


'^) 


^5 


= 2.0 


(for 


>!" ) 


1^6 


= 0.5 


(for 


z) 


0)7 


= 2.0 


(for 


z') 



(6.1) 



The state of convergence of the velocity profile is shown in 
Figure 7(t) ) . 

When a small positive number e is determined, as in Equation 

(5.16), it becomes a parameter for determining the boundary layer 

thickness d. Therefore, attention is focused on how its value is 

selected. The error will be 0.01 if the external end of the boundary 

layer is the point where the physical quantities in the boundary layer 

amount to 99^ of the main jstreaml . . Therefore, a value of '-B^o.oPx 6=0.0006) 

-4 
■'is desirable. Thus, a value of e = 10 will be quite sufficient. 

In actual fact, according to the calculation results, the relation- 

_ _ _4 

ship between £ and d (normalized at d when e = 10 ) will be as shown 



in Figure 8. 'It is only natural that if a large value of|£ is adopt ed,'ithe 



boundary layer will be evaluated with a reduced thickness. Even if 

the value is less than e = 10 , it is known that there will be a 
considerable reduction. In view of the above considerations , -B^«=io-*) 
was used as the criterion for convergence. 

The next question which arises is whether a uniform value can be 
obtained regardless of the departure values (the uniqueness of the 
solution determined by the initial values ) . The departure values /12 
were varied (cj)"^ = 0.0, 0.25, O.5O, 0.75, 1-00; T^ = 0.7, 0.9; z^ = 

0.01, 0.21, 0.41) in order to study the numerical uniqueness. As a 
result of these calculations, it was confirmed that in all cases 



there was a convergence to a uniform solution. Typical examples of 
this are shown in Figure 9. Uniform solutions are obtained even 
though the number of trials differs. It Is possible to select these 
departure values at random, but if the values are selected in the 



25 



1.3 
1.2 

1.1 

1.0 

0.9 

0.8- 

0.7- 

0.6- 



1 1 — I I I I i i | 1 1 — I I I M i | I I I I I I m r I — I I I 1 1 

d/d 

6=10"* 



0.5 

10-* 



III 



I ' ' 1 1 



-I I I Mill 



10- 



10-" 



10" 



10- 



Figure 8. Effects of the convergence criterion 
on the boundary layer thickness 

x=0, 7, = 1000°C, ^»=1.0pcm, ^=latm.| 



vicinity of the anticipated values, naturally this will contribute 
to speeding up the calculations . 



The calculations performed here were for cases when there were 
three unknown boundary conditions (departure values) and three 
constralntjcondltlons . However, It Is known that this method Is valid 
when a similar method Is applied In cases with one or two departure 



values [15] (the work Is much easier ,. since there are no constra int ! 
conditions). Therefore, the general argument In §5.2 Is believed to 
apply also In cases when there are four or more unknown boundary con- 
ditions. (The task of numerical confirmation Is left for future 
study. ) 



As was mentioned In §5.2, It Is desirable to consider the sup- |, 

plementary condition equation (5.8), as this will give more perfect :, 

results as the necessary and sufficient conditions. However, even If ., 

this is omitted, there ought not to be any deficiencies as far as the , 



26 



necessary conditions are 
concerned. Assuming that 
the former is Case 1, and 
the latter is Case 2, 
calculations were per- 
formed under the follow- 
ing conditions : 

1273° K, R^ = 1.00 cm. 



\ = 



p, = 1.0 atm, and the de- 
parture values of (j)" 
^ w 



0.0 



= 0.01, T^ = 0.90 




Departure values 



0.40 



Zic 



0.30 



0.20 



• O^«'=0.00 
A rH.=0.70 

□ Z«r=0.01 

ATh'=0.90 
I Gf 2^=0.41 



0.10 



"W ' w 

When both cases were com- 
pared, the results shown 
in Table 1 (convergence 
values of the departure 
values) were obtained. 
There were no very great 
differences and, as was 
only natural, the calcu- '> 
lating time was less in 
Case 2. The lack of dif- 
ferences in the converg- 
ence values is based on 
the nature of the solu- 
tion, as described in 
§5.2. That is, there are 
no overshoot solutions 
(Figure 5)} and all the 

physical quantities have an asymptotic tendency' m^ving_ towards 0. / 
All the calculations below were performed according to Case 1, taking 
Into consideration the necessary and sufficient conditions. 



4 5 6 

No. of_t rials] 



1 

Figure 9. 

jf=0, Moo=5.74, r,= 1573°K, 5,=0.75cm, p, = {atm. 



Uniqueness of the initial 
values 



/1 3 



A profile of the various physical quantities inside the boundary 
layer obtained in this manner is shown in Figure 10.. If cj) and T, 

■^ WW 



27 



TABLE 1. COMPARISON OF CALCULATION RESULTS 
BY DIFFERENCES IN_THE BOUNDARY CONDITIONS 
AT THE EXTERNAL EDGE) OF THE BOUNDARY LAYER 

[Case I: Including the supplementary condi- 
tions of Equation (5.8)'. Case 2: Not in- 
cluding the supplementary conditions of 
Equation (5-8)] 





Case 1 


Case 2. 


Number 


of trials 


6 




il 5 




5». 




■ -0.02873 




* -0. 02895 


', 


<>'" 




0.79907 




' 0. 79922 




T. 




0.71379 




0.71418 




z. 




0. 05375 




0. 05381 




z'. 




-0.02719 




-0. 02749 




E 




0. 00008 




0. 00007 




V- 




3.7 

132 s ec| 




3. 7 

104 s ec| 


Calculat 
(using 


ing tin 
HITAC 


le, . 1 
5020S)1 



are known, it will be possible to obtain the ablation rate m a nd[ 
the surface temperature T from the following: 



w 



* p.fi.u.r.^ Y' , /, ,2s\ 2s ^ -\ 



(6.2) 



(6.3) 



r.=r,r. I 

Figures|ll(a) and (b) show the effects of the stagnation point tempera- 
ture on the ablation velocity and the surface temperature, plotted 
with the radius of curvature of the object as the parameter. When a 
comparison is made between cases when the thickness of the object 
was made, infinitely great /(y^; = °°) — ^^ thaf'ls, "wTien the object was a 



\ solid hemisphere — and cases when the thickness of the object was J 
flnife "(1/2" of the radius of curvature of the object; y, = 0.5 R, )|^-' •' 

tha-t--is, when the object was a hollow hemisphere with a thickness . - 

half the object's radius of curvature — it is found that. in the " • ^ 
latter case the ablation velocity and the surface temperature are 
rather low, while the heat transfer rate, shown in Figure 11 (d), 
increases on the contrary. This is due to the following fact. 



28 



1.0 

T i' 



0.9 



0.8 



-0.5 — ^ 



0.7 



0.6 



.0.5 -1.0 



■■■ 1 

-0 


1 


1 1 


I 


\ ■* -^ 


/ 


/^ 


V 


\ s^ 


J 






■ \ / 


f 




-0.1- 


1 / I 








\* / r"" 


^_ 






\ / / 


^ 


\X 




\ / / 








\ 1 1 








\l fl 








\ Ik 




^ 


-0.2 - 


\ / /' V 


y 


' - 




-\i.o//' ^ 


<., 





, 


v / 








- yi 






-0.3- 


/A 








// \ / 


_z_ 






//' \ /'^^ 








//' ^i 






. 


III \ 






« "' 


/ \ 






-0.4 - 








1 


\ 






1-2.0 , 


L... 


^"^^sL 


, -0.5 



0.10 
z 



0.05 



Since a uniform tempera- 
ture T, was assumed at 

y = yi_ J in the latter 

case it is necessary to 
maintain a lower tempera- 
ture on the surface of the 
object and to suppress the 
ablation effect. The 
thermal input increases 
for this reason. 

In order to indicate 
the effects of the object 
thickness, calculations 
were made of the ablation 
velocity at the stagna- 
tion point of the object 
and the surface heat 
transfer rate at T, = 

1000° C. The results are 

shown in Figures|12 (a) 

and (b). These results 

indicate that, when the 

object becomes thinner 

at the same thermal input, 

there is an increase in 

the heat transfer rate because of the necessity to suppress the 

ablation velocity; this results in a deterioration of the ablation 

cooling effect. As the object thickness y, becomes smaller, the 

deterioration of the cooling effect proceeds rapidly ahead. Although 
there is little possibility of using an object which is excessively 
thin, it is believed that it would be possible to put into actual'ase 



5 6 
V 



Figure 10., Profile of the physical 
quantities in the boundary layer 

'x=0, ■T,=1500°C, ■S.=1.00cm, ^=iatm.l 



The difference between the surface 

^b 



a hemisphere of R = 1.00 cm. 

heat transfer rate of an object of y, = 0.25 Rj^ '^ 0.5 R^ and that of 



717 



29 



xlO"' 



'(«Mm-''K<!-') 



-1 — I — I — r 



-1 1 1 1 — JT- 

i 



2 - 




=2.00cm 



.J I I I I I ' ' ' 



500 



1000 ISOO 

7.CC) 



Figure 11(a). Ablation rat e 

x=0, Mco=5.74, ^(=latm. 



680 

Twi'C) 

660 



620 



580 



T — I — I — I — I — I — rrr 




_1 I I I I u 



500 



i 1500 



: Figure 11(b). Surface temperature 

*=0, il/oo=5.74, ><=latin.l 



this method (y^ = «') ; 



this method (y, =0.5 



R,); 



[5] (y^ = -) 



0.15 



Z«.- 



0.10 



0.05 



I I I I I 1 1 — I 1 — I r 



y»=0.5;?» 



S»= 0.75cm 




_1_ I I I I I I 1- 



500 



1000 



■r.co 



1500 



5.0 



4.0 - 



3.0 - 



2.0 



1.0 



1 1 1 r 1 I .nfM..,,, 

|7. rn 


1 1 1 ■ 


»» = 0.5R, 


y?.=0.75cm 


/ 


1.00 ^ 


f /" 


^ 


7 


/ 


■■\ — 


/ 


; '^ 


2.00cm 


1/ /'' y "^\^ 


3.00 





500 



1000 



r,ca 



1500 



Figure 11(c) Concentration of Figure 11(d). Surface heat trans- 
ablation gas on surface of object fer rate when there Is a blation 

ir=0, Moo=5.74, /i7=latin.) «=0, Moo=5.74, A=latm.] 



30 



1.1 


1 1 1 1 1 1 1 1 


-1 1 1 


1 1 


1 


- 


8 

,'iio 


- 




^*=oo 


__^-.— . 




O 

<0.9 
■i 
0.8 


"t^ 










V\\/?6= 1.00cm 
W 2.00 












\ 3.00 








0.7 


- ] 








- 


0.6 


- 








- 


0.5 












n 1 


'i 1 " 


1 1 1 


1 1 


1 









l.D 



2.0 



3.0 4.0 _ _ 5.0 
VblRb 



Figure 12 (a). Effects of the thickness of, the 
■ object on the ablation rate • '■'■'- 



^ 1.06 r 

Qo/(Qo) 

1.05 



1.04 
1.03 
1.02 
1.01 
1.00 



*=0, Afoo=5.74, r,= 1273"'K, />,=latni., 



-'-^ 


'■■'■■'■ '■■'■■T-r-1" r 1 — 


1 


' 


- 






- 


-1 


■ 

^(,= 1.00cm 

\ y 2.00 




- 


I // 3.00 


J 


^^ 


• * 


- 


' ' 


^^ 




y;=oo 

\ 


• 1 1 1 


1 


-« J 

1 







1.0 



2.0 



3.0 



4.0 _ ,5 5.0 



Figure 12(b). Effects of the thickness 
of the object on the surface heat trans- 
fer rate 

x=0, M=o=5.74, f;=1273''K, ^=latinj 



31 



0.4 



0.3 



0.2 



0.1 



T — I — I — I — I — I — I — I — I — r 



gt=3.00cni 



T^ 




2.00 



1.00 



7^ 



0.75 




i?»=0.5/?» 

J \ I I I I I L_l I I ' 



500 



1000 1500 

TA°C) 



Figure 13Ja)_; Heat protection 

effectl,of ablationi 



i=0, Moo=5. 74, /'< = latm. 



1000 



(calgf-') 




f.CC) 



Figure 13(b). Cooling effect 
of ablation 

jr=0, M» = 5.74, /!r=latm, J,= oo| 



an object of y, = °° (a solid hemisphere) is only about 6^ at most. 

This is attributable to the good thermal insulation properties of the 
ablation substance Itself. 

Figure 11(c) gives the concentration of the ablation gas on the 
object surface, which was calculated at the same time. It is clear 
that there is the same tendency as with the ablation velocity. 

Figuresil3(a) and (b) show the cooling effect of ablation. First, 
Figure 13(a) shows the ratio between the heat 'quantity ^ ~ ^h 

which is protected by ablation and the heat quantity q when the-re — - 
is no ablation. It is clear that it increases together with the 
increase of the stagnation point temperature of the main fistreamjj]f| and ' 
that the cooling effect improves at higher temperatures. The fact 
that the cooling effect is better in a solid hemisphere than in a 



32 



hollow one points towards a reasonable result, namely, that a larger 
ablation effect gives a better cooling effect. If we define the 
cooling effect of ablation as : 



A=-^^ 



Hw fft9 



(6.4) 



the value of this will display results of the same type as those 
above, as is clear from the broken lines in Figure 13(b). When cal- 
culations are made of the effective heat of ablation, which is often 
used as an index of the ablation effect [17] 3 



H,„=—^ (cal.gr-i) 



ffto 



(6.5) 



this also gives the same results as those for A, as is shown by the 
solid lines in Figure 13(b). 

The fact that there is a superior cooling effect when there is 
a bigger radius of curvature can be explained in terms of the follow- 
ing circumstances. The boundary layer in this case becomes thicker. 



and there is a larger heat quantity to be absorbed] inside the 
boundary layer -i with the^identTcal heat input] ' from the main stream'.] 
The slight differences from the results reported in the literature 
[5] are attributable to the manner of evaluating Equation (4.4). If 
Equation (4.4) is further refined, it is thought that this difference 
will become smaller. The uniqueness of the solutions obtained from 
this direct solution method is guaranteed also by the fact that there 
are no big differences between these solutions. 

Ablation is accelerated when there is an increase in the stagna- /I 9 

tion point pressure of the main current. In Equation (6.2), since 

^ "I 

V is in direct proportion to the - -^ power of p , [refer to Equa- 
tion (A.23)], the following will apply if the other conditions are 

uniform: m«i^» «r*«'^.»^[ , Since p.»«^| , in the final analysis the abla'tToh ^_ 
velocity will be equal to the 1/2 power of the stagnation point 
pressure. 



33 



The first problem arising In numerical calculations by the H-W 
method Is the effect of Ax, the step] amplitude in the x direction. 
This effect Is shown In Figure l4. The ablation velocity at various 
points (x = 0.20, 0.40,. 0.60) was 
calculated for three different step] 
amplitudes: Ax = 0.20, 0.10, and 
0.05. The results were compared 
with the hypothetical ablation velo- 
city at Ax ^ 0, and It was found 
that there was an error of only 
about 0.3^ at the most at Ax = 0.05- 
As can be calculated by Equation 
(3-1). When Ax = 0.05, there will 
be the following value for a' step 



1.1 

thilHa 
1.0 



0.9 - 



amplitude In the n direction of 
n = 0.1: 



0.8 



1 


— 1 1 r—" 

X =0.20 
/ 0.40 


1 


/x =0.60 

1 1 1 







0.05 0.10 



0.15 0.20 
Ax 



0.25 



A./(A,)2<M3^^ 



(6.6) 



Figure l4 
;ude 

7i = 970°C, ^.= 2. 83cm. 



Effects of stepl 
amplitude Ax^ln H-W method 



This also fulfills the stability 
conditions for the difference solu- 
tion of a parabolic partial differential equation and Is believed to 
be suitable with reference to the numerical treatment. Since s = 

1 4 

•jj- X near the stagnation point, it is true that As = 0.003 even when 

2 
Ax = 0.10, and the stability conditions are satisfied when As/(An) 

1 



0.3 < 



2* 



However, the conditions become unstable downstream from 



this. In view of the above results, the calculations were carried 
out with Ax = 0.05. 



In Figures 15(a), (b), and (c), the ratio between the ablation 



rate ml(x) in the flow direction calculated in this way and that at 
the ^stagnation point is shown. The test values and the solutions 
obtained by the microdisturbance theory derived from [5] are also 
shown for the sake of comparison. The present solutions were seen 
to have an even better correspondence with the test values. Stable 
solutions can be obtained only up/ to about x = O.Tj and it is 



/20 



34 



1.2 I 1 1 \ 1 r- 



Mlcrodlstu rba n ce| 
theory I'd^Tl 



0.5 



'This metho_dl 




o Test_ value [5] | 



J I I I L 



0.2 0.4 0.6 0.8 1.0 1.2 

X 



1.2 

1.1 

tk/tko 

1.0 



0.9 - 



0.8 - 



0.7 - 



0.6 



0.5 



1 1 1 1 1 


-^ c -^-5-lJ ■ - 


Thlsll^\N/\ 
•method! /^\ \ 


m=m„cos(fjy^ \ 

\ 


o'Test value [5] . 

1 1 1 r 1 



0.2 0.4 0.6 0.8 1.0 1.2 



Figure 15(a) . Vari ation of] 
of ablation ve'locity 

(a) iWoo = 5.74, f.==970°C,- 1 
Rt=2.S3cm, p,= latm. 



Figure 15(b). Variation of 



V of ablation velocity 

(b)Mc»=5.74, t;=1045°C, 

S"»=1.40cm, ^, = latm, aj;=0.05 



impossible to obtain any solutions 
downstream from this point on ac- 
count of the essential validity of 
the hypothesis and the accumulation 
of numerical errors . The model 
used in these calculations is that 
of a solid hemisphere (y, = <») . 

The hypothesis of semi- 
stationary ablation which was used 
in developing the theory holds that 
ei'^her the amounts of the varia- 
tions do not depend on the time or 
that the time-clependence is ex- 
tremely small. The ablation velo- 
city may justifiably be considered 
a function of x solely. As is 



1.2 

It 

I 1.1 

1.0 



0.6 



0.5 



T 1 1 r 



Thls_method| 
[5] 




m = TWoCos 



o.Test value [5]| 



J -. I I I L 



0.2 0.4 0.6 0.8 1.0 1.2 

X 



Figure 15(c). Variation of 
of ablation velocity 



(c) Moo=5.74, 7-«=1120''C./?. = 1.15cm, 
^, = latm, 6j:=0. 05. 



35 



shown In Figure l6 , It has been observed in experiments [l8] that, 
after the elapse of a definite time, the relationship between the 
surface sweepback distance 7] and the time t becomes more or less 



linear, and that moc-^ 



Is either unrelated to the time or has only 



(2) 

an extremely small dependence. ^- In the case of a perfectly steady 

ablation r|l9], the following would apply: 

m(i-) = m.co9(-|^)j (5^7) /21 

The values are shown graphically for various different cases. The 

slight differences are due. It Is believed, to the dependence on 
time. Within the range of lower stagnation point temperatures In 



the main .stream] , the thermal Input Itself Is found to be low In the 
experiments, while there Is little heating at all points except those 
In the vicinity of the stagnation point In the object Itself. There- 



fore, the distribution of the surf ace__s_we_epb_ackJveloclty becomes 
smaller as one moves downstream from the stagnation point, and It 
becomes more and more difficult for the seml-statlonary hypothesis 
to remain valid. 

To sum up the results of the calculations, the ablation velo- 
city at X = 0.6 can be expressed as In Figure 17 with the radius of 
curvature of the object as the parameter. The distribution In the 
X direction of the heat transfer rate and the surface temperature 
are also calculated In a similar manner. 

The question at issue here is the conduction of heat in the x 
direction Inside the object. Within the range of main current stag- 
I natljDn point temperatures handled here, this has a value which is 
negligibly small in comparison with that in the y direction [5], but 
within the higher temperature ranges it becomes doubtful whether 
this is valid or not.. Popper, Toong, and Sutton [20] have recognized 



^^^The fact that I"! Increases in the negative direction at first 
in Figure 16 indicates that there is elongation due to the thermal 
expansion of the model. 



36 



6.0 



^ 5.0 



4.0 



i 3.0 



2.0 



I.O 



0.0 



-1.0 








■ — r 

Mode 


\ 












y 






..^y<X\ 








V. . 


y. 




<y 


V-^y^/^ 


\ 




:/ 


y 






^ — -^ 


10 


1 1 

150 200 250 
r(sec) 



Figure l6. Measured values of 

surface sweepback distance for 

Teflon [18] 

;Ma.=5.74, f, = latm T,= 1200''C, ^=1.21cm.| 



hlrria 


1 1,1 1 


1 


1.1 


^6= 3.00cm 




1 n 


/ 2.00 


- 




III 0.75 


0.9 


^^^====^r^ 


0.8 


^^^^^^^ 


^^ 


0.7 


^ 


^> — N 

X 


0.6 


II 1 1 


1 



1000 1200 1400 

T.CC) 



Figure 17- Ablation velocity 
at X = 0.6 



A/<»=5. 74, ^ = latm. 



the results obtained with the mlcrodlsturbance theory of the writer 
and colleagues within the low temperature range and have obtained 
good results with tests at higher temperatures using beryllium and 
graphite as the materials. Their tests Included the three-dimensional 
heat conduction effect. 

7. Conclusion 

The direct solution method for sublimating ablation Is one of 
the most effective methods of thermal protection for objects re- 
entering the atmosphere at hypersonic speeds . This method was 
studied In a unified manner from the following three standpoints. 

\ 

(1) Essentially, the ablation phenomenon ought to be deter- 
mined by a combination of the chemical reaction process of the sur- 
face substance and the aerodynamic heat transfer process in the ; 



37 



circumference of the object. The numerical findings obtained here 
indicate that, if the conditions inside the object and the conditions 
of the main streamj are given, it will be possible to solve the 
boundary layer equation by adopting the chemical properties of the 
surface substance as the boundary conditions on the surface of the 
object and also by using the heat conduction inside the object as 
the boundary conditions in the form of the energy balance on the 
surface of the object. By solving the equation in this manner, it 
Is possible to obtain categorically the physical quantities such as- 
the ablation velocity or the surface temperature. 

(2) The Hartree-Womersley method was applied. In this method, 
the differential term of one variable in a parabolic partial differ- 
ential equation is replaced by a difference term, and the equation 
is reduced to an ordinary differential equation. The solution for 
points up to X = 0.7 downstream from the stagnation point was ob- 
tained as in (1) above. Since the distribution of the ablation velo- 
city in the x direction is smooth, the stability conditions are not 
very rigid when the differences are taken. ^ -- •- 

(3) Since the number of boundary conditions on the wall' does 
not satisfy the necessary number requirements, the ordinary differen- 
tial equation obtained in (2) will be a two point boundary value 
problem. One of the boundary conditions on the wall is known. Three 
of the remaining six are completely independent', and the other three 
are related in terms of constraint] conditions. The presence of . 
constraint] conditions makes the problem more difficult. In this re- 
port, the conventional initial value method was further developed, 
and a least square error method was proposed at the edge] of the 
boundary layer. It was found numerically that this method can be 
applied generally without regard for the number of unknown boundary 
conditions or the number of constraint| conditions. The problem of 
sublimating ablation was taken as an example to prove this. 

Calculations were performed by the above-described method, 
using Teflon as the ablation substance. The calculation conditions 
were the following: main istreamj Mach number 5-74, main stream] 



38 



stagnation point temperature 800° C - 1500° C, main ' Stream| stagna- 
tion point pressure 0.5 atm =2.0 atm, radius of curvature of the 
object 0.75 - 3.00 cm, and temperature inside the object of 20° C. 
The following conclusions were reached from the calculation results^ 

(i) Within the range of the : calculations performed here, the 
heat protection effect of ablation improves as the main jstreaml ' 
stagnation point temperature increases. This indicates the advan- 
tages of the ablation cooling method. 

(ii) When there is a greater radius of curvature of the object, 
the heat protection effect increases because the boundary layer be- 
comes thicker. 

(iii) The heat protection effect is better in solid models 
than in hollow models. However, thanks to the good thermal insula- 
tion properties of the ablation substance itself, the difference is 
not large, and there are possibilities that even substances with 
thin walls amounting to about 25^ of the radius of curvature can.-,-, 
be put into actual use. 

(iv) As for the ablation velocity distribution downstream in /22 
the flow direction, there is a rather good coincidence when the main 
'Stream] stagnation point temperature is above 1000° C, in comparison 
with the solutions obtained by the microdisturbance theory and with 
the test values. This indicates that this method can be applied 
effectively to problems of this type. 

It is expected that even better results could be attained if 
the two-dimensional effect of heat conduction inside the object were 
taken into consideration. However, in this case it would be neces- 
sary to use some method other than the H-W method, and the analysis 

at higher temperature ranges will ^be a task to be solved in the i 

future. 



39 



With reference to the solutions of the two point boundary value 
problems, the author wishes to express his heartfelt gratitude to 
director Horlkawa and specialist Watanabe of the Measuring Department 
for enabling the author to study the behavior of the solutions by 
means of an analog computer. 



APPENDIX A 

Boundary Layer Equations 

The following Is the boundary layer equation written In terms 
of the system of coordinates (X, Y) fixed In space for the flow pat- 
tern Illustrated In Figure 1. 

Mass conservation: 



d(pr.) ^ d {pUrc') 4. d (pV>.) ^Q 
dt SX dY 



Momentum conservation: 



3< gx dY 

dX , dY V gy/ 



'' Ri, SY. 



(A.l) 



(A. 2) 
(A. 3) 



Energy conservation 



dt 3X dY 3j 






(A.il) 



40 



Diffusion: 



^^-^^1|-^# 






(A. 5) 



Here, the quantities which have no suffixes 1, 2 are the mean quan- 
tities of a two-component mixed gas. 



Cp=C,iK+Cp2a-'K) 



(A. 6) 



The concentration is supposed to be K, e K. Since K-. + K^ = 1, it 
is true that Kp = 1 - K. Therefore, the analyses given below have 
been given a unified concentration of K. 

Since sweepback occurs on the surface of the object at a velo- 
city of V, J we can use the following transformation to coordinates 



(Xj y) fixed on the surface of the object: 






(A. 7) 



If we use the semi-stationary hypothesis, we obtain: 



dt Sy St 



1^JL^--v,J- 



9y 



(A. 8) 



By means of these, it is possible to transform Equations (A.l) 
through (A. 5) into a simultaneous equation at the system of coordi- 
nates (x, y ) . Furthermore, we can use the non-dlmensionallzed 
quantities expressed by: 



x=x/Rb, y=y/y*, r,=r,/Rt, 
u=u/u*, v=v/v*, p=p/p,i,, 
T=T/T,,H=H/Cp2T,, p='p/T.7u*, 



(A. 9) 



The- following non-dimensional simultaneous partial differential- 
equation can thus be obtained: 



41 



3 (pur,) y djpvr,) 
Bx By 



3u 



Bu 



p„__+,„_^^ 



=0, 
Bx 



3 / Bu \ 
ByK'' By ), 



Sp _ N 



(A. 10) 
(A. 11) 
(A. 12) 



3H , „ BH 

pu— — +pv—-— 

Bx By 



B 
By 



BH 



Pr By 
1 M 3"') 



) 



-{(X-,),o<c„-c„>rf:} 



By 



BK , ,3K S / ^3K\ 



p=fiRT, 
H=C,T+ i Wtt\ 



(A. 13) 
(A.liJ) 

(A. 15) 
(A. 16) 



Let us assume that the shock wave layer is ;lnejllastlc and has a 
constant .density rj Since it is assumed that the velocity component 
on the surface of the object which is obtained from the inelastic 
hypersonic flow theory (Hayes and Probstein [21]_)[wlll be equal to u 

the velocity component at the end of the boundary layer. 



e' 



«,=«» V-y-yrfr''"'* 



(A. 17) 



Here 






(A. 18) 



Let us suppose that 



In this case. 



Furthermore, since 



'"^^^-- 



»,=sin*=a,(x)J 



H.=Cp.T.+ \Wu.^\ 



under conditions where H ' = .|1, 



r.=q^(i--tM'«.2)=r.w 



(A. 19) 



(A. 20) 



jA,_.2i)_ ■-■ 



(A. 22) ' 723 



and the conditions at the boundary layer end will be given, 



42 



It is sufficient for the quantities to be non-dimensionalized 
as follows :| 



y* = (/<.» Rt/p.M*) ' V* =>* u*/St 



(A. 23) 



in order for the original boundary layer equation to be non-dimen- 
s ional and to assume the validity of Equations (A. 10) - (A.l4).| 
Using the state of the mean strean)],, let us rewrite the right side 



of y in the following manner 



,r(2rAfco2-(r-l))^ kr-\) Moo^+i\KR,y (r+i) Mo,-,, 
i^=(;;«/p,»««/?.)'r-!^-— ■ -L_L:_^ Y 



V-^/(i+« 



• J (2rA/co2-(r-l)}^{(r-l)M«2+2}^ /?. 
"^ (r+l)MJ{8r/3(l+-i)2)^ 



(A. 24) 



V R. 



Furthermore , 



ifc, w=u*ycp2T,= y-^-/{i+x)Y 



~<M«?aoc?ICp2T, 



(A. 25) 



(A. 26) 



In this case, we can use a,=ri?r,,c/.2/Cf2=r, c/»2-CD2=i?| to rewrite it as 
follows : 



(A. 27) 



43 



APPENDIX B 

Transfer Coefficient of Two-Component Gas 

Let us suppose that the first approximation of the Chapman- 
Enskog theory [22] applies to the two-component mixed gases which 
are being considered here. Let us then apply the Lennard-Jones (6, 
12) potential to the Inter-partlcle forces. 



(1) Coefficient of viscosity 



/i=/<i/(l+Gi2-;i^ 



Ml 1-A" 



^)+/<2/(l + G2,5J^^J3^)| 



(B.l) 



Here J 






(B.2) 



y-,-^an^ Up are the coefficients of viscosity of the component gases. 



expressed In the following terms: 

/i(=266. 93X 10-' (M,f ) * /ff,2j?,(2. 2>(gr.cm-i.sec-i) 



B.3) 



This Is a relationship which was derived for monoatomic gases, but 
it is known that it can be applied also to multlatomic gases as_w^J:.l_ 
[23]. Here, a. represents the radius of collision (A) and fi . ^^.'__n2__: 

the non-dimensional collision Integral. These are quantities which 
accompany the Lennard-Jones intermolecular function^ [ If we substi- 
tute Equation (B.3) in Equation (B.l) and convert it into non- 
dimensional terms, we obtain: 

+ (l-z)/(l+(mG-l)z}](r/T.)5^ (B.4) 



Here, 



G = G21 , 11= m/fz, m = Mi/ Ml 



(B.5) 



44 



It Is assumed that the component gases present in the same space 
have the same temperature. 



(2) Thermal conduction coefficient 

7: M. K 



Here , 






(B.6) 



(B.7) 



K-, and Kp are the thermal conduction coefficients for the component 



gases , 



/c, = -l^ -^ ft, (cal.cm-i.sec-'.oK-i) 

4 Mi 



(B.:8) 



The Eucken factor Eu. Is a coefficient for approximating the thermal 
conduction coefficients of multlatomlc gases. It Is: 



Eu,=Q. 115+0. 2UCPiM,/R \ ( B . 9 ) 

When these are converted Into inon-dlmenslohal terms, we obtain: 



Here , 



+ (l-z)/{l + (mS-l) ^}](Yr.)' 



(B.IO) 



/(0. 115* +0. 354M2C/>j) j 
5=1.065G2, 



(B.ll) 



45 



(3) Chapman-Rubesln number 



/24 



The Chapman-Rubesln number Is defined as : 

C=pil/p7it. \ (B.12) 

In view of Equation (3-9) and Equation (B.4), It can be expressed 
as follows : 






(b:13) 'I 



(4) Prandtl number 



Pr = ClflKr, Pr, Cp-!^/~ 



(B.14) 



Here, Prp Is the Prandtl number for the air at the external end of 
the boundary layer. In view of Equation (B.l4), 



Pr = Pr2 {].-o.z)\jz /{llG + a-fiG)z] 
+ a-2)/(l + (mG-l)j] 
/[«/(pG+ (l-?5)zj+(l-2) 
/(n-(mG-])2)] 



(B.15) 



(5) Schmidt number 



Sc-'/ifpD 



=Sc.(:-(i-^)^}-£L(4.). 



(B.16) 



Here, SCp Is the Schmidt number of the air at the external end of the 
boundary layer. The two-component diffusion coefficient D Is 



(cm^.sec-') 

Ol2 = C<'l + 02)/2 (A), 



(B.17) 



46 



The non-dlmenslonal collision integrals were used for Q^ ^— 



In view of Equation (B.17), 

(B.18) 



+ a-2)/{n-(mG-l)rj] 



References 

1. Sutton, G. W. The Hydrodynamics and Heat Conduction of a 
Melting Surface. J. Aero/Space Sci., Vol. 25, 1958, pp. 
29-32, 36. 

2.- V Robert sVL-I Stagnation-Point Shielding by Melting and , ' 
Vaporlzation^l NASA TR R-10, 19 59- 

3. Steg, L. and H. Lew. Hypersonic Ablation, The High Temp'erature 

Aspect of Hypersonic Plow. Pergamon Press, 1964, pp. 629-680 

4. Swann, R. T. and J. South. A Theoretical Analysis of Effects 

of Ablation on Heat Transfer to an Arbitrary Axisymmetric 
Body. NASA TND-741, I96I. 

5. Karashima, K., H. Kubota and K. Sato. An Aerodynamic Study of 

Ablation Near the Region of Stagnation Point of Axially 
Symmetric Bodies at Hypersonic Speeds,;] ISAS Rept . No. 425, 
Inst. Space and Aero. Sci., Univ.]Tokyo, I968. 

6. Mills, A. F., A. V. Gomez and G. Strouhal. The Effect of Gas 

Phase Chemical Reactions on Heat Transfer to a Charring 
Ablator-;, AIAA Paper No. 70-869, 1970. 

7. Nomura, Shlgeakl._. Abureshon ni yoru yodomi-ten kuriki kanetsu 

no gensho no'.sgkuitei (Measurement of Reduction of Stagnation 
Point Aerodynamic Heating by Means of Ablation). Koku Uchu 
Gijutsu Kenkyusho Hokoku (Technical Report of National 
Aerospace Laboratory), TR-I67, I968. 



Hartree, D. R. and J. R. Womersley. A Method for the Num'ericalf 
or Mechanical Solution of Certain Types of Partial Differen- 
tial Equations. Proc. Royal Soc . , I6IA, 1937, PP- 353-366. 



~\ ;.l 



^ 47 



9. Smith, A. M. 0. and D. W. Clutter. Machine Calculation of 

Compressible Laminar Boundary Layers.] AIAA J., Vol. 3, No. 4, 
1965, pp. 639-647. 

10. Nachtheim, P. R. and P. Swlgert. Satisfaction of Asymptotic 

Boundary Conditions In Numerical Solution of Systems of 
Nonlinear Equations of Boundary-Layer Type. NASA TND-3004, 
1965. 

11. Fox J L., ed. Numerical Solution of Ordinary and Partial 

Differential Equations. Pergamon Press, 1962. 

12. Lewallen, J. M. A Modified Quaslllnearlzatlon Concept for 

Solving the Nonlinear Two-Point Boundary Value Problem^ NASA 
TND-4025, 1967. 

13. Bartlett, E. P., W. E. Nlcolet and J. T. Howe. Heat-Shield 

Ablation at Superorbltal Reentry Velocities j AIAA Paper 
70-202, 1970. 

14. Rashls, B. and R. N. Hopko. An Analytical Investigation of 

Ablation. NASA TMX-300 , I96O. /25 

15. Kubota, H. An Aerodynamic Study of Hypersonic Heat Shield 

Problem with Local Mass Injection at Multiple Stations. 
ISAS Rept . No. 460, Inst. Space and Aero. Sci . , Univ. Tokyo, 
1971. 

16. Madorsky, S. L. Thermal Degradation of Organic Polymers. 

Polymer Reviews, Vol. 7, Intersclence, 1964, pp. 130-140. 

17. Martin, J. J. Atmospheric Reentry^ An Introduction to Its 

Science and Engineering. Prentice-Hall, 1966. 

18. Karashlma^ Keilchl, Kiyoshl Sato and Hlrotoshl Kubota. 

Kyokucho onsoku jiku-talsho donto buttal no yodoml-ten fukln 
nl okeru abureshon no jlkken (Ablation Experiments in the 
Vicinity of the Stagnation Point in Axlally Symmetrical 
Blunt-Headed Bodies at Hypersonic Speeds). Tokyo Daigaku 
Uchu Koku Kenkyusho Hokoku (Tokyo University, Institute of 
Space and Aeronautic Science, Reports), Vol. 4, No. 3 (A), 
1968, pp. 335-347. 

19. Halns , F. D. Equilibrium Shape of an Ablating Nose in Laminar 

Hypersonic Flow. AIAA J., Vol. 8, No. 7, 1970, pp. 1354- 
1355. 

20. Popper, L. A., T. Y. Toong and G. W. Sutton. Three-Dlmensional 

Ablation Considering Shape Changes and Internal Heat Conduc-- '' 
tlon. AIAA Paper No. 70-199, 1970; also 

Axlsymmetric Ablation with Shape Changes and Internal Heat -> 
ConductionTi AIAA J., Vol. 8, No. 11, 1970, pp. 2071-2074. 



48 



21. Hayes, W. D. and R. F. Probsteln. Hypersonic Flow Theory. 
Academic Press, 196O. 

''22. Dorrance, W. H. Viscous Hypersonic Flow. McGraw Hill, 1962, 

23. Bird, R. B., W. E. Stewart and E. N. Llghtfoot. Transport 
Phenomena. John Wiley & Sons, 19 60. 



Translated for National Aeronautics and Space Administration under 
c'ontract No. NASw 2035, by SCITRAN, P. 0. Box 5^56, Santa Barbara, 

California, 93108 



49